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A nonstandard finite difference scheme 
for solving three-species food chain 
with fractional-order Lotka-Volterra 

model 


S. Zibaei and M. Namjoo* 


Abstract 


In this paper, we introduce fractional-order for a model of tritrophic food 
chain Lotka-Volterra. Moreover, we discuss the stability analysis of fractional 
system. The nonstandard finite difference (NSFD) scheme is implemented to 
study the dynamic behaviors in the fractional-order Lotka-Volterra system. 
Numerical results show that the NSFD approach is easy to implement and 
accurate when applied to fractional-order Lotka-Volterra system. 


Keywords: Fractional differential equations; Lotka-Volterra model; prey- 
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1 Introduction 


Biological systems have been studied for many years. In these systems, it 
is common that state variables represent nonnegative quantities, such as 
concentrations, physical properties, the size of populations and the amount 
of chemical compounds [15]. These biological models are commonly based 
on the systems of ordinary differential equations (ODEs). Exact solutions 
of these systems are rarely in access and usually complicated; hence good 
approximations are required. Numerical methods are often the method of 
choice. They should describe the dynamic behavior of the systems, produce 
the nonnegative solutions, and reproduce the real dynamics of the biological 
systems. The interspecies interaction is among the most intensively explored 
fields of biology. The existance of many mathematical models in that area 
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help us understand the population dynamics of analyzed biological systems. 
Mathematical models of predator-prey systems, characterized by decreasing 
growth rate of one of the interacting populations and increasing growth rate 
of the other, consist of systems of ODEs. In most of the modeled interac- 
tions, all rates of change are assumed to be time independent, which makes 
the corresponding systems autonomous. It is not always possible to find the 
exact solutions of the nonlinear models that have at least two ODEs. It is 
sometimes more useful to find numerical solutions of these types of systems in 
order to programme easily and visualize the results. By applying a numerical 
method on a continuous differential equation system, it becomes a difference 
equation system, i.e., discrete time system. While applying these numerical 
methods, it is necessary that the new difference equation system provide the 
positivity conditions and exhibit the same quantitative behaviours of con- 
tinuous systems such as stability, bifurcation and chaos. It is well known 
that some traditional and explicit schemes such as forward Euler and Runge- 
Kutta are unsuccessful at generating oscillation, bifurcations, chaos and false 
steady states, despite using adaptative step size [13,17,18]. For forward Euler 
method, if the step size h is chosen small enough and the positivity conditions 
are satisfied, the local asymptotic stability for a fixed point is saved while 
in some special cases Hopf bifurcation cannot be seen. Instead of classical 
methods, NSFD schemes can alternatively be used to obtain more qualitative 
results and remove numerical instabilities. These schemes are developed for 
compensating the weaknesses, such as numerical instabilities that may be 
caused by standard finite difference methods. Also, the dynamic consistency 
can be represented by NSFD schemes [10]. The most important advantage 
of this scheme is that by choosing a convenient denominator function instead 
of the step size h, better results can be obtained. If the step size h is chosen 
small enough, the obtained results do not change significantly, but if the step 
size h gets larger this advantage comes into focus. 


As it is well known, in the field of mathematical biology, the traditional 
Lotka-Volterra systems are very important mathematical models which de- 
scribe multispecies population dynamics in a nonautonomous environment. 
Many important and interesting results of the dynamic behaviors for the 
Lotka-Volterra systems have been found in [3, 19,20], such as the existence 
and uniqueness of solutions, the permanence, extinction, global asymptotic 
behavior and bifurcation. Because of the good memory and hereditary prop- 
erties of fractional derivatives, it is often necessary to study the corresponding 
fractional systems. Therefore, the dynamical analysis of the fractional Lotka- 
Volterra systems has attracted a great deal of attention due to its theoretical 
and practical significance. 


Many important results regarding stability of fractional systems have been 
obtained. For instance, the stability, existence, uniqueness and numerical 
solution of the fractional logistic equation are investigated in [7]. The stability 
and solutions of fractional predator-prey and rabies models are discussed 
in [1]. In addition, bifurcation properties of fractional systems have been 
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studied in some papers. For example, conditions for the occurrence of Hopf 
’s bifurcation are explored based on numerical simulations in [29]. The critical 
values of the fractional order are identified for which Hopf ’s bifurcation may 
occur based on the stability analysis in [29]. Thus, it is significant to study 
the dynamical behaviors in the fractional population systems. 


Analysis of fractional Lotka-Volterra equations which are obtained from 
the classical Lotka-Volterra equations in mathematical modeling by the re- 
placing first order derivatives by fractional derivative of order a (0 <a <1) 
have been the focus of recent research in this field. Lots of universal phenom- 
ena can be modeled to a greater degree of accuracy by using the property of 
these evolution equations. The fractional differential equations have gained 
much attention recently due to the fact that fractional order system response 
ultimately converges to the integer order system response. 


The current technological advance has made it possible for humans to dis- 
turb the environmental balance in nature that may cause immense damages, 
such as species extinction or starvation. Therefore, understanding the be- 
haviour of the interaction between the species may help biologists and other 
related parties to prevent those events from happening. The real interaction 
of prey-predator in nature is complex and comprises both interspecies and 
external environmental factors. Therefore, several simplifications are usually 
assumed so that a basic model can be constructed and then developed or 
modified to approach the real system. 


The Lotka-Volterra equations are a system of ODEs in the following form: 
xv’ = ax — bry, 


/ — 


y = -cy+dzy, 
x(0) = 20, y(0) = yo, 


where x and y are prey and predactor, respectively. Here a is the prey 
growth rate in the absence of the predators, b is the capture rate of prey 
per predator, d is the rate at which each predator converts captured prey 
into predator births and c is the constant rate at which death occurs in the 
absence of prey. They show that ditrophic food chains (i.e. prey-predator 
systems) permanently oscillate for any initial conditions if the prey growth 
rate is constant and the predator functional response is linear. 


The classical food chain models with only two trophic levels are shown 
to be insufficient to produce realistic dynamics [5]. Therefore, in this pa- 
per, by modifying the classical Lotka-Volterra model, we analyse and sim- 
ulate the dynamics of a three-species food chain interaction. With non- 
dimensionalisation, the system of three-species food chain can be written 
as 
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x’ = ax — bry, 
y’ = dry — cy — eyz, (1) 
z! = gzy — fz, 
«(0) = 20, y(0) = yo, 2(0) = 20, 
where x, y and z denote the non-dimensional population density of the prey, 
predator and top predator, respectively. The predator y preys on x and the 
predator z preys on y. Furthermore a, b, c, d, e, f and g are the intrinsic 
growth rate of the prey, the death rate of the predator, the death rate of the 
top predator, predation rate of the predator, the conversion rate, predation 
rate of the top predator and the conversion rate, respectively. 

This paper is organized as follows: In the next section, we give some ba- 
sic definitions and properties of the Griinwald-Letnikov (GL) approximation 
and provide a brief overview of the important feature of the procedures for 
constructing NSFD schemes for ODEs. In Section 3, we introduce fractional 
order into the model that describes Lotka-Volterra system and also stability 
theorem and fractional Routh-Hurwitz stability conditions are given for the 
local asymptotic stability of the fractional systems. In Section 4, we will dis- 
cuss the stability analysis of fractional system. In Section 5, we present the 
idea of NSFD scheme for solving the fractional order Lotka-Volterra model. 
Finally in the last section, numerical results show that the NSFD approach 
is easy to be implemented and accurated when applied to fractional-order 
Lotka-Volterra system. 


2 Preliminaries and notations 


In this section, some basic definitions and properties of the fractional calculus 
theory and nonstandard discretization are discussed. 


2.1 Fundamentals of fractional-order 


Fractional differential equations (FDEs) have gained considerable importance 
due to their application in various sciences, such as physics, mechanics, chem- 
istry and engineering [16]. In the recent years, the dynamic behaviors of 
fractional-order differential systems have received increasing attention. Al- 
though the concept of the fractional calculus was discussed in the same time 
interval of integer-order calculus, the complexity and the lack of applications 
postponed its progress till a few decades ago. Recently, most of the dynami- 
cal systems based on the integer-order calculus have been modified into the 
fractional order domain due to the extra degrees of freedom and the flexi- 
bility which can be used to precisely fit the experimental data much better 
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than the integer-order modeling. For example, new fundamentals have been 
investigated in the fractional-order domain for the first time and do not exist 
in the integer-order systems such as those presented in [9, 16]. 


2.2 GL approximation 


The GL method of approximation for the one-dimensional fractional deriva- 
tive is as follows [16]: 


D° x(t) = f(t, x(t)), x(0) = 20, t € (0, ty], (2) 
D° a(t) = jim nv : cap (“Ja (t — gh), 


where 0<a< 1, D® denotes the fractional derivative and h is the step size 
and [4 |] denotes the integer part of a Therefore, Eq. (2) is discretized as 
follows: 


poe =F (ins Ln), S123 ok 


where t, = nh and cH are the GL coefficients defined as: 


l+a 


=i ; 


eras ans j =1,2,3,... 


ay) 


2.3 NSFD discretization 


The initial foundation of NSFD schemes came from the exact finite differ- 
ence schemes. These schemes are well developed by Mickens [13,14] in the 
past decades. These schemes are developed for compensating the weaknesses 
such as numerical instabilities that may be caused by standard finite differ- 
ence methods. Regarding the positivity, boundedness and monotonicity of 
solutions, NSFD schemes have a better performance over the standard finite 
difference schemes, due to flexibility to construct a NSFD scheme that can 
preserve certain properties and structures, which are obeyed by the original 
equations. 


The advantages of NSFD schemes have been shown in many numerical 
applications. Gonzalez-Parra et al. [4] developed NSFD schemes to solve 
population and biological models. Jordan [8] constructed NSFD schemes for 
heat transfer problems. 
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We now give an outline of the critical points which will allow the con- 
struction of NSFD discretizations for ODEs. 


Consider the autonomous ODE given by 


g=f(z), 20)=a,. te [0,é,], 


where f(a) is, in general, a nonlinear function of x. For a discrete-time grid 
with step size, At = h, we replace the independent variable t by 


txt, = nh, n=0,1,2,...,N 
where h = Me The dependent variable x(t) is replaced by 
x(t) © &n, 


where £,, is the approximation of x(t,). 


The first NSFD requirement is that the dependent functions should be 
modeled on the discrete-time computational grid. Particular examples of this 
include the following functions [13, 14]. 


LY © 2n41Yn — Ln41Yn41; 
Os 
LS Ln4+10n,; 


3 


Pare (Ga Pat 


2 a 
A standard way for representing a discrete first-derivative is given by 


Tn+1— In 


h 

However, the NSFD scheme requires that x’ has a more general representation 
Tr4+1— In 
ee 
where the denominator function, i.e. @ has the following properties: 

(i) o(h) =h+O(h?), 

(ii) @(h) is an increasing function of h, 

(iii) }(h) may depend on the parameters appearing in the differential 
equations. 


The paper by Mickens [14] gives a general procedure for determining ¢(h) 
for systems of ODEs. An example of the NSFD discretization process is its 
application to the decay equation 
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where X is a constant. The discretization scheme is as follows [14] 


as Ln Sane, o(h, d) = 1 =e 
Another example is given by 
g! = ye — Az”, 
where the NSFD scheme is 
EE 2 Nias Apmis gee 
) Ay 


It should be noted that the NSFD schemes for these two ODEs are exact in 
the sense that x, = x(t,) for all applicable values of h > 0. In general, for 
an ODE with polynomial terms, 


xv’ =ax+(NL) NL = Nonlinear terms, 


the NSFD discretization for the linear expressions is given by Mickens [14] 


——— =atn +(NL)n, 
where the denominator function is 
eah — | 
h — 
(ha) = 


It follows that if x’ is a function of x which does not have a linear term, then 
the denominator function is just h, i.e. @(h) = h. 


By applying this technique and using the GL discretization method, the 
following relations are yielded: 


n+1 


= 5 Ch in4i—j + f (tnt, En41) 
j=l 
In+1 = oa , n =0,1,2,... 
0 


where cf = o(h)~°. 
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3 Fractional-order Lotka-Volterra model 


Now we introduce fractional-order into the model (1) of Lotka-Volterra 
chaotic system. The new system is described by the following set of frac- 
tional ODEs of order a1, a2,a3 > 0, in the following form 


D“ x(t) = ax — bry, 
D y(t) = dry — cy — eyz, 


D 2(t) = gzy — fz, (3) 
x(0) = 0, y(0) = yo, 2(0) = 20, 
0<a; <1, A= 1,2,3: 


Now, stability theorem on fractional-order systems, fractional Routh- 
Hurwitz stability conditions and their related results are introduced. The 
first stability theorem has been given for incommensurate fractional-order 
systems. 


Theorem 1. ( /12/) Consider the incommensurate fractional-order system 
D°x(t) = f(a(t)), (0) = 20, (4) 


where @ = (Q1,..-,Qn), a; € (0,1] fori = 1,2,...,n and x € R”. The 
equilibrium points of (4), are calculated by solving the equations: 


f(x) =0. 


These points are locally asymptotically stable if all eigenvalues X of the Jaco- 
bian matriz J = gt evaluated at the equilibrium points satisfy: 
xv 
Onn ir 
larg(A)| > se = max(Qy,...,Qn). 


Theorem 2. (HH) Consider the commensurate fractional-order system (4), 
1.€., A, = Q2 . An = aX. If all eigenvalues of the Jacobian matrix of 


an ganar point satisfy: 


arn 
larga) >, 
then, the fractional system is locally asymptotically stable at the equilibrium 
point. 


Consider the system of ODEs given by 


X' = F(X,Y,Z), 
Y' =G(X,Y, 2), (5) 
Z' = H(X,Y,Z), 


A NSFD scheme for solving three-species food chain ... 61 


where F, G and H are nonlinear functions. Let X, Y and Z be the steady- 
state solution, i.e., 


F(X,7,2Z) =G(X,7,2) = H(K,7,Z) =0. 
Now consider small perturbations to the steady-state solutions 


X(t) =X+4+<c(t), 
Y(t)=Y + y(t), 
Z(t) =Z+2(t). 


Frequently these are called perturbations of the steady-state. Substitut- 
ing, we arrive at 


On the left-hand side we expand the derivatives and that by definition 
X= = Z' = 0. 


On the right-hand side we now expand F’, G and H in a Taylor series about 
the point (X,Y, Z). The result is 


a! = F(K,Y,2) + F(X,Y, 20+ B(XY, Dy 
+F,(X,Y,Z)z + terms of order x”, y”, 27, ry, 
yz, xz, and higher, 


y = G(X, Y,2)+ G(X, Y,Z)2+G,(X,¥,Z)y 
+G,(X,Y,Z)z + terms of order x”, y”, 27, xy, 
yz, xz, and higher, 


Z= A(X,Y, Z) + H,(X,Y,Z)x + Hy (X,Y, Z)y 
+H.(X,Y,Z)z + terms of order x”, y”, 27, ry, 
yz, xz, and higher. 

Again by definition, 


F(X,Y,2Z) =G(X,¥ 


s 
ll 
m 
ie 
x 
N 
II 
= 


so we are left with 


/ 

v= a11v a12Y a132, 
oa 

Y = 421L + A22Y + 4232, 
/ 

Zz = a31xr + azgay + 4332, 


where the matrix of coefficients 
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G11 412 413 


A= | aa a22 a23 
431 432 433 12,2 Se 
PLN OY (2): AMY 2) dG Yu) 
ae, LY 2) G,(X,¥,Z) GAX.Y,2Z) |, 
H, ir Y 3) H,(X,Y,Z) H,(X,Y,Z) 


is the Jacobian of the system (5). Hence, the problem has been reduced to 
a linear system, ie., w’ = Aw with w = (x,y, 2)", for states that are in 
proximity to the steady-state (X,Y, Z). 


The Jacobian matrix J of the system (3) at the equilibrium point E = 
(x*,y*, 2") is computed as 


a — by* —ba* 0 
J(E) = dy* c+dx* — ez* ey” d (6) 
0 gz" =F oy. 


The existence and local stability conditions of these equilibrium points are 
as follows: 


Let D(P) denotes the discriminant of a polynomial P 
P()) =X? aA? + aad + ag = 0, (7) 
and 
D(P) = 18a,a2a3 + (a;a2)? — 4a3(a1)° — 4(a2)* — 27(a3)?, 


using the results of [2], we have the following Routh-Hurwitz stability condi- 
tions for FDEs: 


i) If D(P) > 0, then the necessary and sufficient condition for the equilibrium 
point FE to be locally asymptotically stable is a, > 0,a3 > 0,a,a2q — a3 > 0. 


ii) If D(P) < 0, a1 > 0,a2 > 0,a3 > 0, then the equilibrium point FE 
is locally asymptotically stable for a < 2/3. However, if D(P) < 0,a1 < 
0,a2 < 0,a@ > 2/3, then all roots of polynomial (7) satisfy the condition 
arg(X)| < 


iii) If D(P) < 0, a1 > 0,a2 > 0,a1a2 — ag = 0, then the equilibrium 
point F is locally asymptotically stable for all a € [0, 1). 


iv) The necessary condition for the equilibrium point E to be locally asymp- 
totically stable is a3 > 0. 


In the next section, we discuss the asymptotic stability of the equilibrium 
point F of the system (3). 
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A Stability analysis of the model 


To evaluate the equilibrium points of the system (3), let 


ax — bry = 0, 
dxy — cy — eyz = 0, 
gzy — fz=0, 


e 


(5,%,0). All calculations were performed by MAPLE. The local stability 


conditions of these equilibrium points are as follows: 


then the equilibrium points are Eo = (0,0,0), E1 = (0, f,-) and E, = 


(i) The Jacobian matrix (6) at the equilibrium point Ep = (0,0,0) is 


a0 0 
J(0,0,0)=|0-c 0 ], (8) 
00 -f 


with the characteristic equation 
P(A) =)? 441d? + a2 + a3 = 0, 
where 
a,=ft+c-a, az =cf —af —ac, a3 = —fac, 


and D(P) in the above equation is 


Therefore, the eigenvalues of the Jacobian matrix (8) corresponding to the 
equilibrium point Ep are Ay = a, Ag = —c and A3 = —f. 

Clearly, if c 4 f then D(P) > 0. Now, since ag < 0; therefore, based on 
part (i) in Routh-Hurwitz stability conditions, the equilibrium point Ep is 
unstable. 


(ii) The Jacobian (6) at the equilibrium point E, = (0, - —£) is 


ag-bf 9 4g 
g 
Fe fa 0 ef 5 
go o€ g g 
0 - 9 
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where the characteristic equation is 


P(A) = 2 +. a1 A? + a2 + a3 = 0, 


with 
n= il . Be, a2 = —cf, a3 = ae oP, 
and again int as ; 


g 
Here, the corresponding eigenvalues are 


ny = SOE, Glee passer 


Obviously, if g?(a? — cf) + bf (bf — 2ag) # 0 then D(P) > 0. Then as 
a1a2 — a3 = 0; therefore, based on part (i) in Routh-Hurwitz stability condi- 
tions, the equilibrium point FE, is an unstable point. 


(iii) The Jacobian (6) at the equilibrium point E2 = (4, ¢,0) is 


be 
0 marr, 0 
ca ad ea 
J(4. 50) = Db 0 Re ’ (9) 
ag — bf 
0 O 
b 


In this case, the characteristic equation is also 
P()) = M+ ad? tanr+a3 = 0, 


where 


_ _ag— bf (ag — bf)ac 


az = ac, a3 = — 5 : 


and 
4ac(bf (bf — 2ag) + g?a? + b?ca)? 
b+ 
Therefore, the eigenvalues of the Jacobian matrix (9) corresponding to the 
equilibrium point F2 are 


D(P) = 


ms) 
te dy = i fac, MEX Wace. 
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Clearly, if bf (bf — 2ag) + g2a? + b?ca 4 0 then D(P) < 0. Now if bf > ag 
then a, > 0,a2 > 0,a,a2 — a3 = 0 and based on part (iii) in Routh-Hurwitz 
stability conditions the equilibrium point EF, is locally asymptotically stable 
for all a € [0,1). 


5 NSFD for fractional-order Lotka-Volterra model 


For system (3) and applying Mickens scheme by replacing the step size h 
by a function ¢(h) and using the GL discretization method, the following 
equations are obtained: 


n+1 


Oey — 
Ci En41—j = ALn — bLn41Yn, 


2 Ch Ynt1—j = —CYnt1 + dtn4iYn — €Yn+12n; (10) 


y Cp enp1—j = —fenti + G2nYn41- 


Comparing equations (10) with system (3), we note the following: 
1. The linear and nonlinear terms on the right-hand side of the first equation 
in system (3) are in the forms 


LY Ln, —xry % —Ln41Yn- 


2. The linear and nonlinear terms on the right-hand side of the second 
equation in (3) are 


YS —Yn4+1) LY + Ln41Yns YZ —Yn412n- 
3. The linear and nonlinear terms on the right-hand side of the third equation 
in (3) are 
2 —2n41, 2Y & ZnYn+t- 


Doing some algebraic manipulations to equations (10) yields the following 
relations: 
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n+1 
—S oe tn41-j + arn 
x ye ais 
oie co) + byn ; 
n+1 
eats + din41Yn (11) 
_ _J=1 
oor 6, Pees, : 
n+1 
—S Fen 41-5 + 92nYn+1 
at eee 
n+1— ee iF ’ 
where 
co’ = b1(h)™, co” = ga(h)~°?, co* = d3(h)°, 
with [21] 
eh 1 ech — 1 efh 1 
oi(h) = : bo(h) = $3(h) = 
a c f 


Proposition 1. The numerical solutions obtained from system (11) for case 
0<a, <1, t1=1,2,3 satisfy 
Ln > 0 En4+1 > 0 
Yn > O> Yn41 > 0 (12) 
Zn > 0 Zn41 >0 


for all the relevant values of n. 


Proof. Since cj" > 0 and by recursive relation 


Leo j =1,2,3.... 


we have c“' <0, j > 0. Now system (11) shows that relations (12) is 
established. For case a; = 1, 7 = 1,2,3 we should consider the following 
system: 


Tn4+1— in 
4 = ALn — bXn41Yn; 
a 
Yn4t1 — Un _ Ld 
b = CYn+1 741 Tn+1Un — €Yn4+1%n> 
2 
en41 — %n = i 
o = —f2n4i + 92nYn41- 
3 


By solving this system for %n41,Yn+1 and z,41 we conclude that relation 
(12) holds. 
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6 Numerical results 


Analytical studies always remain incomplete without numerical verification 
of the results. In this section, we present numerical simulation to illustrate 
the results obtained in the previous sections. The numerical experiments 
are designed to show the dynamical behaviour of the system in three main 
different sets of parameters and initial conditions: 

(i) The case where bf = ag, 


(ii) The case where bf > ag, 


(iii) The case where bf < ag. 


To show the dynamics of the system (3), set the parameter a = b=c= 
d=e= ff =1 given as fixed parameters and g as a varied parameter. 


(i) The case where bf = ag 

For the case bf = ag the equilibrium point E2 has three eigenvalues with zero 
real part corresponding with stable centre point in xy plane . We consider 
the case a, = @2 = a3 = 1 which corresponds to the classical Lotka-Volterra 
system. Figures 1 and 2 represents the phase portrait for solutions where 
parameter g = 1 with the initial conditions (a(0), y(0), z(0)) = (0.5, 1,2), 
for simulation time 40s and step size h = 0.1 and h = 0.5. In this case, 
prey x, predator y and top predator z persist and have populations that vary 
periodically over time in a common period. 


Once again an equilibrium is achieved within the system, such that each 
predator population increases as the population of its respective prey in- 
creases. Each predator population also peaks and then begins to decrease 
shortly after its respective prey population peaks and begins to decrease. 
The plots of populations x and y are essentially the same as they were in 
the 2D system, and the new predator population z behaves similarly with 
respect to y as y behaves with respect to x. All three populations share a 
common period. 
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x(t), y(t), Z(t) 


Tat 
SM MS 
wu. ov ys 
oN 

rid le Pi t 


‘ Vi ' 
wy ay 
Soh eet a sf al N 
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t 


Figure 1: Plot of populations 2, y and z over time for the case bf = ag with ay = a2 = 
a3 =landh=0.1. 


x(t), y(t), 2(t) 


ie 


it x(t) y(t) 


Figure 2: Plot of populations z, y and z over time for the case bf = ag with a1 = a2 = 
a3 =landh=0.5. 


Figures 3 and 4 depict the phase trajectory of the fractional-order Lotka- 
Volterra chaotic system (3) for commensurate order a; = a2 = a3 = 0.90 and 
parameters g = 1 with the initial conditions ((0), y(0), z(0)) = (0.5, 1, 2), for 
simulation time 40s and step size h = 0.1 and h = 0.5. 


x(t), y(t), Z(t) 
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co 


L 1 1 L L 1 L 
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Figure 3: Plot of populations 2, y and z over time for the case bf = ag with ay = a2 = 
a3 = 0.90 and h = 0.1. 
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Figure 4: Plot of populations z, y and z over time for the case bf = ag with ay = a2 = 
a3 = 0.90 and h = 0.5. 


Figures 5 and 6 depict the phase trajectory of the fractional-order Lotka- 
Volterra chaotic system for incommensurate order and parameters g = 1 with 
the initial conditions (x(0), y(0), z(0)) = (0.5,1,2), for simulation time 40s 
and step size h = 0.1 and h = 0.5. 
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Figure 5: Plot of populations z, y and z over time for the case bf = ag with a, = 


0.99, ag = 0.95,a3 = 0.90 and h = 0.1. 
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Figure 6: Plot of populations z, y and z over time for the case bf = ag with a, = 


0.95, ag = 0.90, a3 = 0.80 and h = 0.5. 


In Figure 7, the phase trajectory of the fractional-order Lotka-Volterra 
chaotic system is depicted for incommensurate order and parameters a = 


1b = 2c =5,d = 4,e = 3,f 


3,9 


6 with the initial conditions 


(x(0), y(0), z(0)) = (0.5, 1,2), for simulation time 40s and step size h = 0.1. 
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Figure 7: Plot of populations z, y and z over time for the case bf = ag with a, = 
0.99, a2 = 0.95,a3 = 0.90 and h = 0.1. 


(ii) The case where bf > ag 

For the case where bf > ag, two eigenvalues for E2 are pure imaginary 
initially-spiral stability corresponding with centre manifold in xy plane and 
one negative real eigenvalue corresponding with stable one-dimensional in- 
variant curve in z axis. Hence, the equilibrium point E>» is locally stable 
spiral sink. On the other hand, prey x and predator y persist and has popu- 
lations that vary periodically over time with a common period. The solutions 
are plotted in Figures 8 and 9 for commensurate order a, = a2 = a3 = 1 and 
parameters g = 0.88 with the initial conditions (#(0), y(0), z(0)) = (0.5, 1, 2), 
for simulation time 100s and step size h = 0.1 and h = 0.5. 
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Figure 8: Plot of populations z, y and z over time for the case bf > ag with ay = a2 = 
a3 =landh=0.1. 
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Figure 9: Plot of populations 2, y and z over time for the case bf > ag with ay = a2 = 
a3 =1landh=0.5. 


Figures 10 and 11 depict the phase trajectory of the fractional-order 
Lotka-Volterra chaotic system (3) for incommensurate order a; = 0.90, 
Q2 = 0.80, ag = 0.70 and parameters g = 0.88 with the initial conditions 
(x(0), y(0), z(0)) = (0.5, 1,2), for simulation time 100s and step size h = 0.1 
and h = 0.5. 
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Figure 10: Plot of populations x, y and z over time for the case bf > ag with ay = 
0.90, a2 = 0.80,a3 = 0.70 and h = 0.1. 
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Figure 11: Plot of populations x, y and z over time for the case bf > ag with ay = 
0.90, a2 = 0.80,a3 = 0.70 and h = 0.5. 


(iii) The case where bf < ag 

For the case where bf < ag, two eigenvalues for HE, is pure imaginary 
initially-spiral stability corresponding with centre manifold in xy plane and 
one positive real eigenvalue corresponding to unstable one-dimensional in- 
variant curve in z axes. Hence the equilibrium point £2 is a locally unstable 
spiral source. In this case, the prey x and top predator z can survive, growing 
periodically unstable. On the other hand, predator y persists and has popu- 
lations that vary periodically stable. The solutions for this case are shown in 
Figure 12 for commensurate order a, = a2 = a3 = 1 and parameters g = 1.6 
with the initial conditions (a(0), y(0), z(0)) = (0.5, 1,2), for simulation time 
50s and step size h = 0.1. 
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Figure 12: Plot of populations 2, y and z over time for the case bf < ag with ay = a2 = 
a3 =landh=0.1. 


In Figure 13 the phase trajectory of the fractional-order Lotka-Volterra 
chaotic system (3) is depicted for commensurate order ay = a2 = a3 = 
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0.50 and parameters g = 1.6 with the initial conditions (a(0), y(0), z(0)) = 
(0.5, 1,2), for simulation time 50s and step size h = 0.5. 
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Figure 13: Plot of populations 2, y and z over time for the case bf < ag with ay = a2 = 
a3 = 0.50 and h = 0.5. 


The solutions for this case are shown in Figure 14 for incommensurate 
order and parameters g = 1.6 with the initial conditions (x(0), y(0), z(0)) = 
(0.5, 1,2), for simulation time 50s and step size h = 0.1. 
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Figure 14: Plot of populations x, y and z over time for the case bf < ag with ay = 
0.60, ag = 0.50, a3 = 0.40 and h = 0.1. 


In Figure 15 the phase trajectory of the fractional-order Lotka-Volterra 
chaotic system is depicted for incommensurate order and parameters g = 1.6 
with the initial conditions (a(0), y(0), z(0)) = (0.5, 1,2), for simulation time 
50s and step size h = 0.5. 
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Figure 15: Plot of populations z, y and z over time for the case bf < ag with ay = 


0.8, a2 = 0.6,a3 = 0.5 and h = 0.5. 


In Table 1 for different step size h, the qualitative results, obtained by 
NSFD scheme, of the fixed point Ez are respectively compared to classical 
methods such as forward Euler and 4th order Runge-Kutta. From Table 1, 
it follows that the CPU time of the method NSFD is less than the CPU time 
of the forward Euler and Runge-Kutta methods. Also if step size h is chosen 
small enough, the results of the proposed NSFD scheme are similar with the 
results of the other two numerical methods. But if the step size h is chosen 
larger, the efficiency of NSFD scheme is clearly seen. 


Table 1: Qualitative results of the equilibrium point FE for different time step sizes, t= 
0-200 for the case where bf = ag 


h Euler CPU time Runge-Kutta CPU time NSFD CPU time 
0.001 Convergence 0.016342 Convergence 0.032029 Convergence 0.000206 
0.01 Convergence 0.014760 Convergence 0.028096 Convergence 0.000205 
O.1 Convergence 0.013917 Convergence 0.027959 Convergence 0.000203 
0.2 Divergence - Convergence 0.025959 Convergence 0.000202 
2 Divergence - Divergence = Convergence 0.000201 
10 Divergence - Divergence - Convergence 0.000201 


In Figure 16 the numerical solution of forward Euler and fourth order 


Runge-Kutta methods are compared with NSFD scheme graphically. 
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Figure 16: Numerical solutions for forward Euler and fourth order Runge-Kutta and 
NSFD methods with h = 0.1 for the case bf = ag and aj = ag = a3 = 1. 


7 Conclusion 


In this paper, we study the fractional-order Lotka-Volterra model. The sta- 
bility of equilibrium points is studied. Numerical solutions of these models 
are given. The reason for considering a fractional order system instead of its 
integer order counterpart is that fractional order differential equations are 
generalizations of integer order differential equations. Also using fractional 
order differential equations can help us to reduce the errors arising from the 
neglected parameters in modelling real life phenomena. 


We argue that the fractional order models are at least as good as integer 
order ones in modeling biological, economic and social systems (generally 
complex adaptive systems) where memory effects are important. 
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